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Q ' Abstract 

Based on a regularized Volterra equation, two different approaches for numerical differ- 
entiation are considered. The first approach consists of solving a regularized Volterra 
_^ \ equation while the second approach is based on solving a disretized version of the 

regularized Volterra equation. Numerical experiments show that these methods are 
efficient and compete favorably with the variational regularization method for stable 
_^ calculating the derivatives of noisy functions, 

j^ | Keywords: ill-posed problems, numerical differentiation. 
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m ■ 1 Introduction 

> 

Calculating the derivatives of noisy functions is of prime importance in many applications. 
The problem consists of calculating stably the derivative of a smooth function / given its 
noisy data fs, \\fs — f\\ < <5. This is an ill-posed problem: a small error in / may lead 
to a large error in /'. Many methods have been introduced in the literature. A review 
is given in [7]. Divided differences method with h = h(S) has been first proposed in [4], 
see also [HI El E] • Necessary and sufficient conditions for the existence of a method for 
stable differentiation of noisy data are given in [8l chapter 15], see also [9]. In our paper 
a method for stable differentiation based on solving the regularized Volterra equation 



Au(x) + f s (0) := / u(s)ds + fs(0) = fs(x), (1) 

Jo 

is proposed (see also [1QJ H [9] ) • One often applies the Variational Regularization (VR) 
method 

\\Au - f s \\ 2 + a\\u\\ 2 ->min (2) 

for stable differentiation. 

In this paper (and in [lj) an approach, based on the fact that the quadratic form of 
the operator A is nonnegative in real Hilbert space L 2 (0, a), a = const > 0, is used. 
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2 Methods 

Consider two different approaches to solving equation dU). The first approach consists of 
solving directly regularized equation ([1]). The second approach is based on the Dynamical 
Systems method (DSM) and an iterative scheme from [3]. 

2.1 First method 

In pQ , the derivatives of a noisy function fs are obtained by solving the equation 

au a ,s + Au a) s = fs- (3) 

If a = a(5) > is continuous on [0, 5q), d~o > and 

lim a(8) = 0, lim — — = 0, (4) 

5^o 5^o a(S) 

then the following result holds (see [T]): 
Theorem 1 Assume ([!]). Then 

lim \\ux — u\\ = 0, 

where u$ solves ([3]) with a = a(S). 

The solution of ([3|) is: 

fs(x) 



I x f x s 

us(x) = ^exp( )/ ex.p(-)f s (s)ds + 

or a Jq a 



(5) 



a 

This formula and an a priori choice a (5) = 6 /c, where k G (0,1), c is a constant, 
yield a scheme for stable differentiation. When a(8) is known, the problem is reduced 
to calculating integral ([5]). There are many methods for calculating accurately and fast 
integral (|5]) (see e.g. |2j). However, there is no known algorithm for choosing k,c which 
are optimal in some sense. The advantage of our approach is that the CPU time for the 
method is very small compared with the VR and DSM, see Section 13.11 Moreover, one 
can calculate the solution analytically when the function f$ is simple by using tables of 
integrals or MAPLE. 

2.2 An iterative scheme of DSM for solving discretizations of the regu- 
larized Volterra equation 

Another approach to stable differentiation is to use the DSM (see [8]). The DSM yields a 
stable solution of the equation: 

F(u) =Au-f = 0, ueH, (6) 



where H is a Hilbert space and A is a linear operator in H which is not necessarily bounded 
but closed and densely defined. The DSM to solve ([6]) is of the form: 



v! = -u + {T + a(s))~ i A*f, u(0) = u , (7) 

where T := A* A and a(t) > is a nonincreasing function such that a(t) — > as t — ► oo. 
The unique solution to (J7|) is given by 



(8) 



u(t) = n e" 4 + e~ l [ e s (T + a(s))" 1 A* f ds . 
Jo 

An iterative scheme for computing u(t) in (J5J) is proposed in [3]: 

u n+ i = e- /l "u n + (l-e- /l ")(T + a n )~ 1 ^*/«5, h» = Wi-^ 

With ao satisfying 

(5 < \\Auao - Ml < 25, (9) 

one chooses a n and h n as follows: 



1 + t ' 



where 1 < q < 2, to = 0- To increase the speed of computing we recommend choosing 
q = 2. At each iteration one checks if 

0.95 < \\Au n - f s \\ < 1.0015. (10) 

This is a stopping criterion of discrepancy principle type (see [3]). If t n is the first time 
such that (fTUj) is satisfied, then one stops and takes u n as the solution to ©. The choice 
of ao satisfying ([9]) is done by iterations as follows: 

1. As an initial guess for ao one takes ao = i||A|| 2 <5 re z, where 5 re i = tt4it. 

J ' ' lull 

II Aw — fx\\ 

2. If ^ = c > 3, then one takes a\ := 2 (c-D as *^ e nex ^ S uess an< ^ checks if 

condition (|10p is satisfied. If 2 < c < 3 then one takes a\ := ao/3. 

3. If ^ — — = c < 1, then ai := 3ao is used as the next guess. 



4. After ao is updated, one checks if (|T0|) is satisfied. If (|T0|) is not satisfied, one repeats 
steps 2 and 3 until one finds ao satisfying condition (fT0|) . 



Algorithms for choosing ao and computing u n are detailed in algorithms 1 and 2 in [3]. 

3 Numerical experiments 

Numerical experiments are carried out in MATLAB in double-precision arithmetic. In 
all experiments, by u(t), Um(t), mdsm(*) an d u VR(t) we denote the exact derivative, the 
derivatives computed by the first, the DSM and the VR methods, respectively. In this 
section by n we denote the number of points used to discretize the interval [0, 1]. 



3.1 Computing the first derivatives of a noisy function 

Let us compute the derivatives of the function f(t) = sin(7rt) contaminated by the noise 
function e(t) = <5cos(107r£). The derivative of /(£) is f'(t) = ircos(-Kt). To solve this 
problem we use three methods: the first method, based on computing integral ([5]), the VR 
method, and the DSM method, based on a discretized version of (pQ). Numerical results 
for this problem are presented in Figured! In our experiments, since the results otained 
by the DSM and the VR are nearly the same, we present only the results for the DSM in 
Figure Q] and [2] in order to make these figures simple. 

In this experiment the trapezoidal quadrature rule is applied to integral equation ([1]) 
and is used for computing integral ([5]). One may use higher order intepolation methods 
to compute integral ([5]). However, it does not necessarily bring improvements in accuracy. 
This is so because using a high order intepolation method for inaccurate data may even 
lead to worse results. This is the case when the noise level is large. 

The approximate derivative formula ([5]) for t close to does not use much information 
about fs- Thus, we only use © for computing f'(t) for t G [f,l]- For t G [0, |)], we 
take gs(t) := fs(l — t) and use formula © for g$(t) with t G (i, 1]. That is, we have a 
discontinuity at t = « of the solution, obtained by the first method in Figured] and EJ The 
same idea is applied in discretizing equation ^) in the implementation of the DSM and 
VR. 

In the DSM and VR we also use the trapezoidal quadrature rule to discretize equation 
(HJ). Since the right-hand side fs contains noise, using high order collocation methods 
does not necessarily improve the accuracy. Experiments have shown that the use of higher 
order collocation methods leads to linear algebraic systems with larger condition numbers 
and yields numerical solutions with low accuracy. 
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Figure 1: Numerical results for f$(x) = sin(7rt) + <5cos(107rf). Discretization points n 
100. 



The CPU times for the VR and DSM are about 0.0125 sec. The CPU time for the 
first method is much smaller: 0.0015 sec. Here, we should bear in mind that the DSM 



and the VR use iterations to look for "good" regularization parameter a while the code 
based on the first method does nothing to look for a but uses a as an input value. If one 
also uses the regularization parameter as an input in the VR and DSM, although these 
methods still take more time than the first method the difference in computation time is 
not so large. 
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Figure 2: Numerical results for fs(x) 
n = 100. 



sin(27rt — ^tt) + <5cos(107rt). Discretization points 



The error of the first method for 5 = 0.02 is larger than those of the VR and the 
DSM, but when 5 = 0.002 then the first method gives smaller errors. From Figure [Hand 
[2j one can see that the solutions obtained by the DSM are better than those obtained 
by the first method for all t € [0, 1] except for the t which are close to the boundary of 
the interval. Indeed, it can be showed analytically that the solution u to equation ([2]) 
satisfies n(0) = u(l) = 0. However, the derivative of / in Figure Q] satifies /'(0) = n and 
/'(l) = — vr. If the computed derivatives at the points close to the boundary are discarded, 
then in both cases the DSM and the VR are more accurate than the first method. 

Figure [2] presents the numerical experiment for f(t) = sin(27rt — ^tt) contaminated by 
the same noise function e(t) = e)cos(107r£). For this problem, since the function to be 
differentiated / satisfies /'(0) = /'(l) = both the DSM and the VR give more accurate 
results than the first method. 

From Figure Q] and [2] one can see that for 5 = 0.02 the computed derivatives are 
very close to the exact derivative at all points except for those close to the boundary in 
Figure [H 



3.2 Computing the second derivatives of a noisy function 

Let us give numerical results for computing the second derivatives of noisy functions. The 
problem is reduced to an integral equation of the first kind. A linear algebraic system is 



obtained by a discretization of the integral equation whose kernel K is Green's function 

K(8,t): 



s{t 
t(s 



if 
if 



s < t 

s > t 



Here s, t G [0, 1] and as the right-hand side / and the corresponding solution u one chooses 
one of the following (see [3] ) : 

s — s 
case 1, f(s) = — - — , u(s) = s, < s < 1, 



case 2, /(s) 



6 ' 

sin(27rs) 

4vr 2 



+ s — 1, u(s) = sin(27rs), < s < 1. 



Collocation method is used for discretization. This discretization can be improved by 
other methods but we do not go into detail. We use n = 10, 20, ..., 100, and b n ^ = b n + e n , 
where e n is a vector containing random entries, normally distributed with mean 0, variance 
1, and scaled so that ||e n || = <5 re j||6 n ||. This linear algebraic system is mildly ill-posed: the 
condition number of j4ioo is 1-2158 x 10 4 . 

Table 1: Results for case 1 and 2 with 5 re i = 0.01, n = 20,40, ..., 100. 
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0.0479 


100 


5 0.2956 
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0.2948 


100 


4 0.0254 


6 


0.0379 



Table Q] shows that numerical results obtained by the DSM are more accurate than 
those by the VR. Figure [3] plots the numerical solutions for these cases. The computation 
time of the DSM in these cases is about the same as or less than that of the VR. From 
Table Q] one can see that both the DSM and the VR perform better in case 2 than in 
case 1. Note that the regularized equation to solve for second derivatives in this case is 
of the same form as equation ([2]). As we discussed earlier, it is because in case 2 we have 
/'(0) = f(l) = 0. 

We conclude that in this experiment the DSM competes favorably with the VR. 

Looking at Figure [3] case 1, one can see that the computed values at t = and t = 1 
are zeros. Again, the regularized scheme forces the computed derivative u to satisfy the 
relations u(l) = n(0) = 0. If one wants to compute the derivative of a noisy function on 
an interval by the proposed method, one should collect data on a larger interval and use 
this method to calculate the derivative at the points which are not close to the boundary. 



4 Concluding remarks 

In this paper two approaches to stable differentiation of noisy functions are discussed. The 
advantage of the first approach is that it contains neither matrix inversion nor solving 
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Figure 3: Plots of solutions obtained by DSM, VR when n = 100, 5 re i = 0.02. 

of linear algebraic systems. Its computation time is very small. The drawback of the 
method is that there is no known a posteriori choice of a(<5). The second approach is an 
implementation of the DSM. It competes favorably with the VR in both computation time 
and accuracy. The DSM competes favorably with the VR in solving linear ill-conditioned 
algebraic systems. A posteriori choice of a, an efficient way to compute integral ([5]) for 
the first method, and an efficient discretization of the Volterra equation (JTJ) with the 
implementation of the DSM are planned for future research. 
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